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INTRODUCTION 


The  objective  of  the  research  described  in  this 
report  was  to  predict  the  spectrum  of  Rayleigh  waves 
generated  by  a  nuclear  explosion  over  an  island  atoll, 
as  recorded  at  teleseismic  distances.  Our  approach  was 
a  computational  one,  since  this  problem  cannot  at 
present  be  handled  theoretically:  generation  of  the 
seismic  waves  takes  place  within  a  small  distance  of 
ground  zero  (which  is  assumed  to  be  located  near  the 
middle  of  the  atoll) ,  and  this  implies  that  the  coupling 
of  energy  from  the  explosion-generated  shock  wave  in 
the  atmosphere  into  seismic  waves  in  the  earth  is 
characteristic  of  the  shallow  water  depth  typical  of 
the  interior  of  atolls  of  the  type  we  are  considering. 

On  the  other  hand,  a  major  influence  in  shaping  the 
spectrum  of  teleseismic  Rayleigh  waves  is  the  deep 
ocean  structure  which  surrounds  the  island,  through 
which  the  Rayleigh  waves  propagate  most  of  the  way  to 
the  recording  stations.  The  transition  from  shallow  to 
deep  water  takes  place  well  within  the  critical  distance 
for  Rayleigh  wave  generation,  so  even  the  approximate 
coupling  techniques  developed  by  McGarr  and  Alsop  (1967) 
and  McGarr  (1969)  are  inapplicable.  The  basic  difficulty 
is  that  at  the  present  time  theory  cannot  handle  laterally 
inhomogeneous  structures  with  the  necessary  precision,  and 
it  is  obviously  wrong  to  assume  a  laterally  homogeneous 
oceanic  structure  between  the  island  source  and  the 
continental  station. 

Consequently,  we  thought  of  applying  numerical 
techniques  to  carry  the  shock  wave  from  the  atmosphere 
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through  the  shallow  water  and  the  island  structure  into 
the  deep  ocean:  the  distributions  of  elastic  displace¬ 
ments  and  stresses  with  depth  completely  specify  the 
input  to  a  well-known  problem  which  has  a  unique  solu¬ 
tion,  viz.,  Rayleigh  wave  propagation  in  a  laterally 
homogeneous  horizontally  plane- layered  medium. 

Since  the  seismic  stations  are  located  on  continents, 
we  anticipated  being  able  to  use  McGarr's  theory  for 
predicting  what  happens  at  the  continental  margin.  An 
alternative,  once  programs  for  numerical  calculations 
of  elastic  wave  propagation  in  an  inhomogeneous  media 
were  developed,  would  be  to  model  numerically  the  Ravleigh 
wave  propagation  through  the  ocean/continent  boundary. 

We  anticipated  that  this  would  be  a  fairly  straight¬ 
forward  project.  Delineation  of  the  problem  is  easy, 
material  properties  of  the  type  of  structure  in  question 
are  known,  it  is  quite  simple  to  set  out  a  procedure 
for  attacking  the  problem,  and  after  all,  numerical 
calculation  of  elastic  wave  propagation  in  solids  has 
been  intensively  studied  for  years;  for  example  by 
Holzer  and  his  colleagues  (1966),  and  by  Alterman  and 
her  colleagues  (Alterman  et  al.,  1970). 

We  were  wrong,  in  that  it  now  appears  that  the 
problem  is  beyond  the  state  of  the  art  of  numerical 
analysis,  and  our  project  ended  before  we  were  able 
to  advance  that  particular  frontier  sufficiently  to 
solve  the  atoll  problem. 

We  started  with  the  most  elementary  approach:  an 
explicit  finite  difference  scheme  (in  both  time  and 
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space)  in  which  the  derivatives  in  the  wave  equation 
were  approximated  by  finite  difference  expressions.  As 
a  starting  point  we  used  the  program  descrioed  by 
Alterman  in  a  number  of  papers,  which  Professor  \lterman 
generously  made  available  to  us,  but  it  soon  became 
clear  that  that  program  could  not  easily  be  converted 
to  the  cylindrical  symmetry  we  required. 

We  therefore  programmed  a  finite  difference  code 
in  which  cylindrical  symmetry  was  built  in.  The 
inhomogeneous  elastic  wave  equation  is  rather  compli¬ 
cated,  and  for  checking  out  the  analysis  we  found  it 
very  useful  to  have  access  to  a  FORMAC  program, 

(Tobey ,  1967).  This  list  processing  system  was  able 
to  take  in  the  multi-dimensional  tensor  equations  and 
produce  the  final  difference  approximation  in  FORTRAN 
language,  with  which  we  could  compare  and  verify  our 
own  deviation  and  programming. 

It  is  prudent  to  approach  computers  with  the  same 
apprehension  and  care  that  bullfighters  use  in  their 
work,  so  we  started  with  simple  test  cases  whose 
solutions  are  known  analytically.  It  soon  became 
apparent  that  there  are  inherent  difficulties  in  the 
finite  difference  method,  namely: 

1.  The  inhomogeneities  of  material  properties  in 
the  problem  are  abrupt:  air  to  water,  water  to  rock, 
crust  to  mantle;  in  addition,  these  material  inhomo¬ 
geneities  are  two-dimensional  (i.e.,  r  and  z  variation) 
and  clso  the  source  function  is  discontinuous  in  time. 
These  factors  appear  to  cause  unavoidable  instabilities 
in  the  numerical  calculations. 
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2.  The  source  function  is  applied  to  one  of  the 
boundaries  on  which  material  properties  are  discontinuous. 
This  also  appears  to  cause  instability. 

3.  Even  though  we  used  a  high-speed  computer  with 
ten  tape  units,  32K  words  of  core,  ah'  a  large  disc, 

it  was  very  difficult  to  arrange  to  make  the  mesh  large 
enough  to  avoid  spurious  reflections  from  the  artificial 
boundaries  caused  by  terminating  the  mesh:  the  P  waves 
went  out  and  back  before  the  Rayleigh  waves  ever  got 
to  the  deep  ocean. 

C.  Constantino  and  F.C.  Karal  discussed  these 
questions  with  us  at  length  (private  communication, 

1970).  The  concensus  was  that  the  difficulties  are 
indeed  inherent  in  finite  difference  approximations  to 
dymanic  problems.  We  subsequently  found  that  even  an 
implicit  scheme,  rather  than  an  explicit  one,  did  not 
solve  these  difficulties. 

We  therefore  decided  to  use  a  finite  element  scheme 
rather  than  a  finite  difference  scheme.  The  finite 
element  method  has  been  used  for  many  years  in  static 
problems  of  structural  mechanics.  In  this  method  the 
material  is  divided  into  small  homogeneous  polygonal 
elements  which  are  welded  together  alcng  their  edges. 

The  elastic  wave  equation  is  solved  analytically  in 
ec.ch  individual  element  under  the  single  assumption 
that  the  displacement  varies  linearly  along  the  edges 
of  the  elements.  Basically,  the  method  works  because 
the  displacements  at  the  vertices  of  each  element,  as 
calculated  in  erch  of  the  contiguous  elements,  must  all 
be  the  same.  This  fact  leads  to  condition  equations 
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which  are  easily  solved  if  a  large  enough  computer  is 
available . 

Unfortunately,  even  after  very  considerable  effort 
we  were  unable  to  prodv.ce  correct  answers  for  elementary 
problems  in  which  shock  waves  are  applied  to  a  material 
discontinuity.  A  detailed  discussion  is  given  in  the 
following  sections.  Further  work  on  this  problem  is 
currently  a  part  of  a  Ph.D.  program  in  the  Department 
of  Geophysics  at  the  Pennsylvania  State  University. 


FINITE  DIFFERENCE  APPROACH 


The  finite  difference  approximation  to  numerical 
solution  of  elastodynamic  problems  is  well  known;  see, 
for  example,  references  given  by  Alterman  et  al.  (1970). 
In  this  approach  the  derivatives  cf  the  wave  function 
with  respect  to  time,  and  the  spatial  derivatives  of 
the  material  properties,  are  approximated  by  finite 
differences . 

In  this  section  we  discuss  our  development  of 
this  approach  in  some  detail. 

The  elastic  wave  equation  is: 


2 

(X  +  y)  V  ;V*S)  +  (7* S)  (VX)  +  y  t  S 

+  +  V(Vy»S")]  =  pQS> 

Here  the  displacement  jT  *  (u,v,w)  and  X,  y  and  p  , 

the  elastic  constants  and  density  respectively,  are 
functions  of  position.  Neglecting  second  derivatives 
of  the  elastic  constants  one  obtains  the  following  for 
the  u  component  equation  in  rectangular  coordinates: 
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where  X  =  X(x,y,z),  u  =  u(x,y,z),  pQ  «  PQ(x,y,z)  and 
the  displacement  in  the  ex  direction.  The  equations  of 
motion  for  the  v  and  w  displacement  component,  (e  and  e 
direction.0 ) ,  have  the  same  form  as  equation  (1) . 

These  equations  in  cylindrical  coordinates  are: 
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where  X  =  X(r,<J>,z),  p  =  p(r,<J>,z),  pQ  =  P0(r,*,z)  and 

u,v  and  w  are  the  displacements  ir  the  e  ,  en  and  e 
, .  .  r’  G  z 

directions  respectively . 

The  cases  we  are  presently  interested  in  are  those 
having  axial  symmetry,  in  which  case  the  equations  reduce 
to: 
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(3u  +  3wj  +  3ji  2.9w 
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We  notice  that  the  displacements  in  the  e  and  e 
...  r  z 

directions  are  completely  uncoupled  to  displacements  in 

the  eQ  direction. 

For  an  explicit  finite  difference  solution  of  these 
equations  one  approximates  the  derivatives  by  either 
backward,  centered,  or  forward  differences  (except  for 
the  derivatives  of  the  elastic  constants,  which  will  be 
discussed  later).  In  regions  other  than  the  free 
surface  (to  be  discussed  later)  it  is  generally 
accepted  that  centered  differences  usually  give  the 
most  accurate  results. 

Thus,  with  r  -  mAr,  z  -  nAz  and  t  *  pAt,  where 
m,n,  and  p  are  integers,  the  derivative  approximation 
would  be  of  the  form: 


r 
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Substitutions  into  equations  (2)  leads  to  an  explicit 
solution  of  the  form: 
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and  a  similar  form  for  the  w  component.  Thus,  knowledge 
of  displacements  at  times  pAt  and  (p-l)At  in  a  neighborhood 
of  locations  about  the  point  (m,n)  enables  one  to  calcu¬ 
late  displacements  at  time  (p+l)At. 

Stability  condition 

The  selection  of  increments  Ar,  Az,  and  At  is  not 
arbitrary,  since  the  difference  equations  must  satisfy 
a  stability  condition,  which  the  case  of  the  wave 
equation  insures  that  information  does  not  flow  faster 
in  the  finite  difference  scheme  than  it  does  in  the 
medium  being  modelled.  Stability  conditions  can  be 
derived  for  the  given  form  of  the  wave  equation  using 
a  number  of  different  methods,  all  of  which  involve  a 
great  deal  of  calculation  and  ultimately  only  give 
approximate  stability  conditions.  We  use  instead  a  new 
method  which  gives  the  exact  stability  condition  quite 
easily  for  the  elastic  wave  equation:  setting  the 
coefficients  of  the  u(m,n,p)  term  equal  to  zero  yields, 
in  the  appropriate  limits,  the  approximate  stability 
conditions  obtained  by  the  other  methods. 

For  the  axisymmetric  case  the  approximate  stability 
condition  is: 


At  < 


-13- 


where 


■('  hi 


If  we  look  at  the  finite  difference  expression  for  the 
u  component,  where 
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We  see  that  setting  the  coefficient  of  the  u(m,n,p)  term 
equal  to  zero  yields  the  stability  condition: 


2  2 

at  <  [(— )“  (i  *  — 1 — j)  ♦  (— ) 
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which  reduces  to  the  earlier  expression  when  AZ  =  Ar  and 
m-*-®. 

For  the  w  component  of  displacement  one  has  a 
different  stability  condition.  The  finite  difference 
expression  for  the  w  component  is: 
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which  gives  for  the  stability  condition 
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These  relations  were  checked  by  running  sample 
cases,  and  the  additional  terms  included  found  to  be 
valid,  i.e.,  they  do  set  the  limit  on  the  point  at 
which  instabilities  develop.  Changes  in  the  elastic 
constants  enter  these  expressions,  and  we  see  that 
knowing  the  exact  stability  condition  is  very  important 
when  dealing  with  abrupt  velocity  changes,  since  at 
certain  boundaries  or  corners  non-stable  effects  may 
develop  and  gradually  contaminate  the  rest  of  the 
problem. 

Method  for  handling  velocity  changes  without  specifying 

boundary  conditions 

If  one  uses  the  homogeneous  form  of  the  elastic 
wave  equation  to  solve  a  problem  in  which  there  are 
several  velocity  discontinuities,  solutions  are 
obtained  by  matching  boundary  conditions  for  each  point 
at  which  there  is  a  velocity  change.  This  approach  is 
impractical  for  complicated  velocity  distributions.  If 
one  uses  the  heterogeneous  elastic  wave  equation  (as  we 
do  in  this  report)  the  elastic  parameters  and  their 
space  derivatives  are  input  and  the  boundary  conditions 
are  automatically  satisfied  within  the  equation. 

One  might  expect  that  the  spatial  derivative  of 
X  and  m  could  be  handled  in  the  same  manner  as  for  the 
derivatives  of  displacement  (using  central  differences); 
however,  this  method  gives  incorrect  results.  A  method 
which  does  work  involves  spreading  velocity  changes 
over  two  grid  points  with  the  boundary  point  having  a 
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velocity  (v^  +  where  and  are  the  velocities 

of  the  adjoining  velocity  zones.  One  then  takes  the 
derivatives  of  the  elastic  parameters  X  and  y  to  be  zero 
at  the  mesh  points  before  and  after  the  boundary  point 
and  approximates  the  derivative  at  the  boundary  by  a 
centered  difference.  This  procedure  is  illustrated  in 
Figures  2  for  the  case  y  =  0. 

If  this  procedure  is  used  while  holding  the 
stability  factor  and  time  increment  At  constant  (or 
the  stability  factor  and  Ar  or  Az,  according  to  the 
direction  of  the  velocity  change),  excellent  trans¬ 
mission  and  reflection  amplitudes  are  obtained.  Trans¬ 
mission  and  reflection  results  are  shown  in  Figure  3 
for  the  problem  in  which  the  P  velocity  varies  in  one 
dimension  from  1  to  2,  and  y  =  0.  In  general  both 
shear  and  longitudinal  waves  are  present.  One  now  needs 
to  pick  X  and  y  values  at  the  boundary  such  that  both 
the  longitudinal  and  shear  velocities  make  a  linear 
transition  from  one  medium  to  the  other.  We  first  pick 
yB,  the  value  of  y  at  the  boundary,  to  be 


yB 
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(assuming  pQ  to  be  unchanged),  thus  setting  the  shear 
velocity  at  the  boundary.  Then  one  can  set  the  longitudinal 
velocity  at  the  boundary  by  setting  the  value  of  X  at  the 
boundary : 

1  _  _  2 
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This  can  be  extended  to  the  case  where  p,  X,  and  y  all 
change  at  the  boundary.  First  we  pick  a  linear  change  in 
p: 


pB 


pl*p2 
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Then  we  pick  a  linear  B  change 


pB  = 
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Thus  we  can  satisfy  the  condition 
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Thus  one-dimensional  velocity  changes  in  a  medium  can 
be  handled  by: 

a)  spreading  velocity  changes  over  two  mesh  points, 

b)  picking  values  of  p,X  and  y  at  the  boundary 
point  such  that  the  velocities  make  linear  transitions, 
and 

c)  setting  the  derivative  of  the  elastic  parameters 
equal  to  zero  except  at  the  boundary  point  and  using  a 
central  difference  there. 

This  method  is  very  effective  for  cases  in  which 
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there  is  inhomogeneity  in  only  one  direction,  but 
unfortunately  cannot  be  extended  to  cases  with  two-  or 
three-dimensional  velocity  variations  (private  communi¬ 
cation,  C.  Costantino). 

Free  surface 


Accounting  for  the  free  surface  is  another  aspect 
of  the  problem  in  which  the  centered  difference  approxi 
mation  is  not  satisfactory.  Consider  the  mesh  defined 
in  Figure  4. 

To  take  the  derivatives 
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one  would  expect  to  be  able  to  write 
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However,  Alterman  and  Rotenberg  (1969)  showed  that 
this  method  becomes  progressively  more  inaccurate  after 
several  time  steps.  They  determined  that  one  should  use 
a  centered  difference  for  the  r  derivative,  but  a  back¬ 
ward  difference  for  the  z  derivative: 


and  the  same  fos  the  w  component  of  displacement. 

u(r,0,t)  and  w(r,0,t),  the  displacements  at  the 
ficticious  layer,  are  easily  evaluated  when  the 
boundary  stress  condition  is  applied.  For  axisymmetric 
problems  it  is 
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The  difference  expression  for  these  conditions  are 
then,  with  a  driving  pressure  P(M,l,t), 


Ar 


from  which  u(r,0,t)  and  w(r,0,t)  can  be  eval  .ated  and 
substituted  into  the  expression  for 


az 


z=0 


and 


which  are  then  used  in  the  homogeneous  elastic  wave 
equation  for  the  mesh  near  the  f^ee  surface.  Thus,  in 
the  case  of  the  free  surface  we  treat  the  boundary 
condition  separately,  using  the  homogeneous  equation, 
whereas  the  other  velocity  changes  are  handled  auto* 
matically  using  the  heterogeneous  equation. 


Non-reflecting  boundary 

Elimination  from  the  observed  signal  of  energy 
reflected  from  the  boundary  introduced  into  the  problem 
by  truncating  the  medium  is  accomplished  either  by 
placing  the  boundary  sufficiently  far  from  each  point 
of  observation,  or  by  devising  a  non-reflecting 
boundary.  Expending  the  mesh  requires  more  computer 
storage  and  longer  run  times,  which  makes  this  alterna¬ 
tive  a  poor  choice. 
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A  non- re  fluting  boundary  in  one  dimension  is  easy 
to  create  because  one  knows  the  relation  between  the  wave 
striking  the  boundary  and  its  reflection,  and  it  is  there¬ 
fore  possible  to  cancel  the  reflected  wave  at  future  times 
near  the  boundary.  This  is  even  possible  to  do  with  both 
P  and  S  waves  arriving,  since  one  can  work  with  both 
transverse  and  longitudinal  displacements.  In  two  and 
three  dimensions  this  simple  approach  becomes  impossible 
and  only  very  poor  approximate  cancelling  schemes  can  be 
devised  in  special  cases. 

Another  possible  approach  would  be  to  introduce 
damping  into  the  equations.  By  increasing  the  damping 
coefficient  while  adjusting  the  elastic  parameters  (to 
keep  the  impedance  of  the  damped  and  non-damped  regions 
equal)  it  should  be  possible  to  remove  all  reflections. 

.  The  easiest  way  to  introduce  damping  is  by  adding 
k'u  and  k'w  velocity-dependent  terms  to  the  u  and  w 
component  equations,  respectively.  This  method  works 
fine  for  a  few  time  steps  and  then  errors  start  to 
accumulate,  causing  the  damping  to  become  inefficient 
and  impractical.  A  more  effective  way  of  introducing 
damping  is  to  treat  the  elastic  parameters  as  operators 
and  to  add  a  time  deiivative 
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One  sees  that  in  the  one-dimensional  case  this 
amounts  to  adding  a  velocity-dependent  term  giving  the 
familiar  damped  response.  In  the  three-dimensional 
axisymmetric  problems  this  approach  would  generate  very 
complicated  expressions  for  the  damped  heterogeneous 
elastic  wave  equation  but  it  supposedly  leads  to  stable 
damping  of  the  wave  train. 
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FINITE  ELEMENT  METHOD 


We  can  summarize  the  drawbacks  of  the  finite 
difference  method  as  follows: 

1.  The  lack  of  a  good  procedure  for  creating  a 
non-reflecting  boundary  where  the  mesh  is  terminated 
makes  it  necessary  to  use  a  very  large  mesh  (the  reflec¬ 
tions  travel  faster  than  the  Rayleigh  wave  we  are 
interested  in)  or  elre  to  rescale  the  material  zones  as 
time  advances. 

2.  The  discontinuities  in  the  atoll  problem  are 
two-dimensional,  and  this  situation  cannot  be  handled 
satisfactorily.  Explicit  codes  are  intrisically  unstable 
for  this  type  of  problem. 

We  therefore  decided  to  use  the  finite  element 
method:  the  evolution  of  this  approach  is  described  in 
the  Introduction  of  the  report. 

Application  of  the  finite  element  method  to  dynamic 
problems  is  reasonably  straightforward.  In  roughly  the 
same  manner  as  the  implicit-explicit  finite  difference 
schemes,  the  static  equations  are  solved  at  each  time 
step,  and  the  stresses  and  displacements  throughout  the 
whole  mesh  are  then  input  as  initial  conditions  to  the 
identical  static  problem  at  the  next  time  step. 

It  turns  out  that  there  are  many  annoying  procedural 
problems  in  creating  an  efficient  program  for  carrying 
out  the  dynamic  finite  element  scheme,  the  most  important 
of  which  is  contriving  to  enumerate  the  vertices  in  the 
mesh  so  that  the  coefficient  matrix  --  which  must  be 
inverted  in  the  program  at  each  time  step  --  is  a  band 
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matrix  of  small  bandwidth.  It  is  not  practical  to  invert 
a  1,000  by  1,000  matrix  at  each  of  1,000  time  steps,  at 
least  until  the  millenium  of  ultra-parallel  computations 
arrives . 

We  made  a  literature  search,  using  services  of  DDC, 
and  found  a  highly  sophisticated  finite  element  program 
--  the  SLAM  code  --  which  was  developed  for  the  Air 
Force  by  the  Illinois  Institute  of  Technology 
(Costantino,  1968).  The  program  documentation  consists 
of  five  inch-thick  volumes,  and  there  was  a  time  delay 
while  our  project  personnel  mastered  them.  However,  the 
total  delay  was  clearly  much  less  than  would  have  been 
required  to  develop  an  equivalent  program  from  scratch. 

Additional  difficulties  were  encountered  at  this 
point.  When  after  a  nine-month  delay  following  our  request 
a  copy  of  the  program  arrived,  we  found  that  it  was  an 
old  version  written  for  the  CDC  6600  and  it  was  riddled 
with  errors,  both  elementary  FORTRAN  errors  and  logical 
mistakes.  It  was  clear  that  what  we  received  was  not  a 
working  version  of  the  code,  and  we  were  unable  to 
obtain  a  working  version  despite  repeated  requests. 

Part  of  the  program  was  written  in  6600  assembly 
language,  so  the  process  of  correcting  the  program 
required  day-to-day  use  of  that  type  of  machine.  We 
obtained  an  in-house  terminal  to  a  6600,  and  in 
addition  we  simultaneously  simplified  and  condensed  the 
code  into  a  version  suitable  for  running  test  uses  on 
the  CDC  1604  which  at  that  time  was  available  in  the 
Seismic  Data  Laboratory  (SDL)  computer  center. 
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After  study  and  experimentation,  ve  decided  to 
adopt  a  generalization  of  the  implicity/explicit  scheme 
used  in  the  SLAM  code,  the  Newmark  Beta  Method  (Newmark, 
1959) .  This  scheme  has  the  advantage  that  it  allows  the 
user  to  control  the  way  the  time-step  acceleration 
varies  between  the  time  steps.  The  method  is  described 
by  three  difference  equations: 
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where  0  is  the  parameter  which  controls  the  variation 
of  acceleration,  M  is  the  mass  mafix,  D  the  damping 
matrix,  and  K  the  stiffness  matrix.  In  the  SLAM  code 
as  we  received  it,  8  ■  0.  This  scheme  is  explicit  only 
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when  D  -  0,  since  in  that  case  successive  values  of 
displacement  and  acceleration  can  be  obtained  by 
inverting  only  the  diagonal  mass  matrix.  The  scheme  is 
explicit  in  general  if  the  velocity  damping  matrix  is 
diagonal. 

The  first  major  revision  we  made  to  the  code  was 
to  install  the  general  implicit  scheme,  which  is 
represented  by  the  "solved"  difference  equations: 
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Numerical  solution  of  these  equations  involved  inverting 
the  augmented  stiffness  matrix  at  every  time  step.  This 
was  accomplished  through  standard  inversion  techniques; 
the  matrix  was  initially  block  upper- triangularized  and 
stored  on  tape.  The  back  substitution  process  was  then 
carried  out  from  this  tape  at  every  time  step  (for 
general  discussion  of  these  techniques,  see  Wilkison, 
1965).  The  storage  required  for  the  matrix  calculations 
necessarily  reduced  the  maximum  mesh  bandwidth  from 
.lOO  to  45  on  a  machine  with  32K  words  of  core.  The  same 
procedure  was  also  programmed  for  a  360/67,  and  the 
advantage  in  computing  speed  over  the  1604  compensated 
for  the  increased  number  of  calculations  required  at 
each  time  step. 

At  that  time,  we  also  considerably  streamlined  the 
SLAM  program  by  deleting  the  option  for  handling 
plastic  deformation,  whj.ch  is  not  relevant  to  our 
problem. 

After  the  program  had  been  checked  out,  we  proceeded 
to  apply  it  to  elementary  problems  whose  solutions  are 
known  analytically. 

The  simplest  of  such  problems  is  an  end- loaded  rod 
constrained  by  rollers  to  have  elastic  displacement  only 
along  its  length.  The  nodes  at  one  end  were  loaded  by 
a  pulse  of  stress  with  a  duration  equal  to  the  time 
step  interval.  The  node  point  displacements  should  thus 
be  step  functions  beginning  at  the  appropriate  elastic 
wave  arrival  time.  The  actual  calculations  showed  oscil¬ 
lations  superimposed  on  the  solution.  Figure  5  shows  the 
calculation.  Replacing  the  puJse  source  with  a  Gaussian 
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time  variation  for  the  source  reduces  these  oscillations 
(Figure  6) .  Introducing  damping  does  not  decrease  the 
amplitude  of  the  oscillations,  but  does  decrease  the 
sharpness  of  compression  in  the  wave  front,  and  makes 
the  back  side  of  the  step  function  into  a  slowly 
decaying  ramp. 

Changing  the  Newmark  B-parameter  or  the  integration 
time  interval  did  not  affect  the  amplitude  of  the 
spurious  oscillations,  but  did  decrease  their  period. 

When  the  rod  thickness  was  doubled  and  divided  into 
two  parallel  rows  of  elements  of  half  the  original 
thickness,  the  calculated  results  were  identical  to 
the  original  rod  with  a  single  row  of  squaTe  elements. 
This  implies  that  slight  changes  in  the  small  main 
diagonal  frequencies,  while  maintaining  the  same  time 
step,  do  not  affect  the  calculations. 

Other  test  cases  involved  varying  the  size  of  the 
elements  in  the  direction  transverse  to  the  axis  of  the 
rod;  none  of  these  variations  affected  the  solutions. 

What  did  affect  the  calculations  was  variation  in  the 
size  of  elements  in  the  direction  of  the  rod  axis. 

Figure  7  shows  the  nodal  displacements  along  the 
rod.  The  amplitude  of  the  spurious  oscillations  increases 
when  the  wave  front  reaches  the  area  of  the  rod  where 
the  element  lengths  change. 

The  wavefront  velocity  is  within  2%  of  the 
theoretical  value,  which  we  regard  as  satisfactory 
agreement. 

We  also  studied  the  classical  problem  of  a  line 
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source  at  the  surface  of  a  homogeneous  half space  (Lamb, 
1964).  The  mesh  used  is  shown  in  Figure  8.  For 
convenience,  and  without  lack  of  generality  for  the 
purpose  of  studying  the  agreement  between  the  calculated 
solution  and  the  theoretically  known  solution,  we 
assumed  a  Poisson  solid;  o  =  1/4.  The  time  step 
was  0.4  seconds,  and  the  P  wave  velocity  1/4  feet/second. 
The  mesh  size  was  chosen  such  that  the  P  wave  travel 
time  across  the  smallest  element  in  the  mesh  was  equal 
to  the  time  step.  A  delta  function  source  was  assumed 
to  be  applied  during  the  first  time  step,  and  removed 
for  all  subsequent  times.  No  damping  was  used. 

The  output  is  shown  in  Figures  9  and  in  the 
frequency-wavenumber  spectrum  (Figure  10) ,  The 
frequency -wavenumber  spectrum  was  computed  from  the 
output  along  the  surface  of  the  halfspace  over  that 
portion  of  the  mesh  where  the  spacing  of  the  nodes  was 
uniform,  i.e.,  0  to  4  feet  away  from  the  source  node. 

The  P  wave  and  Rayleigh  arrivals  are  accurately 
calculated,  but  the  Rayleigh  wave  is  followed  by 
spurious  oscillations.  The  Rayleigh  wave  has  the  proper 
horizontal  and  vertical  amplitudes  and  polarity,  and 
has  retrograde  elliptical  particle  motion. 

In  checking  through  possible  reasons  for  the  errors, 
.ve  discovered  that  the  boundary  condition  at  the  free 
surface  was  not  being  properly  met.  Even  though  the 
forces  at  all  the  nodes  at  the  surface  of  the  mesh  were 
constrained  to  vanish,  the  vertical  and  radially 
tangential  stresses  did  not  vanish. 
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We  therefore  decided  to  constrain  the  boundary 
stress  through  the  use  of  the  stress  tables  already 
calculated  by  the  program.  This  was  accomplished 
through  a  Lagrange  multiplier  technique.  To  illustrate 
it,  we  first  consider  the  static  problem  which  is 
solved  by  ninimizing  the  strain  energy: 

u  =  ix^Kx  (5) 


with  the  stress  matrix  constraint: 


(6) 


The  solution  to  the  problem  posed  by  equations  (5)  and 
(6)  is  obtained  by  inserting  a  vector  of  Lagrange 
multipliers,  X  and  minimizing  the  modified  strain  energy: 
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The  minimization  of  equation  (7)  gives  the  solution  as: 
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same 


where  f  is  now  the  effective  force. 

The  dynamic  problem  proceeds  in  much  the 
manner,  resulting  in  the  following  equations: 
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The  third  version  of  SLAM  solves  these  equations.  The 
penalty  paid  is  an  increase  in  run  time  by  a  factor  of 
two  and  a  decrease  in  maximum  mesh  bandwidth  from  45 
to  29.  This  version  is  still  capable  of  being  run  in  a 
32K  machine.  With  it,  the  two  stress  components  t  and 

Ld  X 

a  are  as  specified  on  the  constrained  surface  for  all 
z  z 

time. 

The  output  from  this  version  of  the  program  is 
shown  in  Figures  11  and  12.  Figure  11  shows  spurious 
oscillations  between  the  P  wave  arrival  and  the  Rayleigh 
wave,  in  addition  to  the  same  "himalayan"  effect  follow¬ 
ing  the  Rayleigh  wave  seen  in  Figure  9. 

The  frequency-wavenumber  spectra  show  that  the 
ringing  is  due  to  travelling  energy  having  lower  group 
velocity  than  the  Rayleigh  wave,  and  which  is  normally 
dispersed.  The  result  of  the  dispersion  appears  to  be 
due  to  the  design  of  the  mesh.  There  also  appears  to 
be  wavenumber  aliasing  in  the  horizontal  direction; 
this  is  also  unexplained.  The  displacement  at  the 
source  node  (Figure  13)  shows  that  these  effects  are 
not  due  to  ringing  at  the  origin.  We  therefore  believe 
that  this  sort  of  mesh  is  not  capable  of  accurate 
representation  of  the  solution,  and  that  more  detail 
in  the  mesh  (possibly  around  the  region  of  the  source) 
would  be  necessary  to  obtain  accurate  results.  However, 
this  refinement  is  beyond  the  capability  of  the  computers 
available  to  us  at  the  present  time. 

One  of  us  (DWM)  has  continued  to  study  this  problem, 
and  has  made  use  of  Wilson's  code  (Wilson,  1969).  He 
obtained  results  similar  to  those  reported  here.  He  has 
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also  modified  the  original  SLAM  code  by  replacing  the 
lumped-mass  matrix  by  a  consistent-mass  matrix 
(Zienkiewicz ,  1967).  The  7esults  of  this  modification 
are  not  yet  available.  The  work  of  Goudreau  and 
Taylor  (1972)  has  provided  new  insight  which  may  aid 
in  solving  this  class  of  elastic  wave  problems. 
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Figure  1.  Schematic  diagram  of  the  finite  difference 
displacement  calculation. 
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Figure  5.  Transmission  and  reflection  of  a  pulse  at  a  material 
boundary  using  the  finite  difference  method. 
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Figure  5.  Longitudinal  displacement  of  the  "standard"  one 
dimensional  rod  with  square  wave  pulse  input. 
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Figure  6.  Longitudinal  displacement  of  the  "standard 
dimensional  rod  with  Gaussian  shaped  pulse  input. 
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Figure  8.  Homogeneous  isotropic  halfspace  mesh  used  fov 
solving  Lamb's  problem. 


96.95- 


Figure  9b.  w  displacement  for  Lamb's  problem  with  unconstrained 
surface  stress. 
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Figure  11a,  u  displacement  for  Lamb's  problem  with  constrained 
surface  stress. 


